function [part] = partialLagvar_LogGMM(beta2, bpartial, mean_a, mean_b, mean_la, mean_lb, mean_ha, mean_hb, U, RAND, n, k)
global draws

% This file estimates the partial effects at the mean (zero for binary variables)
% Tobias Pfutze July 2009


% Compute matrix of drawn errors ER
for j=1:draws
    PART(1:k-2,j)=(bpartial(1:k-2)*exp(mean_a*beta2(1:k)+U(RAND(j),1))...
        +(bpartial(1:k-2)-bpartial(k:2*k-3))*exp(mean_a*beta2(1:k)+mean_b*beta2(k+1:2*k)+U(RAND(j),1)+U(RAND(j),2)))...
        /((1+exp(mean_a*beta2(1:k)+U(RAND(j),1))+exp(mean_b*beta2(k+1:2*k)+U(RAND(j),2)))^2);
    
    PART(k:2*k-3,j)=(bpartial(k:2*k-3)*exp(mean_b*beta2(k+1:2*k)+U(RAND(j),2))...
        +(bpartial(k:2*k-3)-bpartial(1:k-2))*exp(mean_a*beta2(1:k)+mean_b*beta2(k+1:2*k)+U(RAND(j),1)+U(RAND(j),2)))...
        /((1+exp(mean_a*beta2(1:k)+U(RAND(j),1))+exp(mean_b*beta2(k+1:2*k)+U(RAND(j),2)))^2);
    
    PART(k-1,j)=(bpartial(k-1)*exp(mean_a*beta2(1:k)+U(RAND(j),1))...
        +(bpartial(k-1)-bpartial(2*k-2))*exp(mean_a*beta2(1:k)+mean_b*beta2(k+1:2*k)+U(RAND(j),1)+U(RAND(j),2)))...
        /((1+exp(mean_a*beta2(1:k)+U(RAND(j),1))+exp(mean_b*beta2(k+1:2*k)+U(RAND(j),2)))^2);
    
    PART(2*k-2,j)=(bpartial(2*k-2)*exp(mean_b*beta2(k+1:2*k)+U(RAND(j),2))...
        +(bpartial(2*k-2)-bpartial(k-1))*exp(mean_a*beta2(1:k)+mean_b*beta2(k+1:2*k)+U(RAND(j),1)+U(RAND(j),2)))...
        /((1+exp(mean_a*beta2(1:k)+U(RAND(j),1))+exp(mean_b*beta2(k+1:2*k)+U(RAND(j),2)))^2);
    
    %Replace discrete changes for binary variables
    
    PART(3:5,j)=exp(mean_ha*beta2(1:k)+U(RAND(j),1))...
        ./(1+exp(mean_ha*beta2(1:k)+U(RAND(j),1))+exp(mean_hb*beta2(k+1:2*k)+U(RAND(j),2)))...
        -exp(mean_la*beta2(1:k)+U(RAND(j),1))...
        ./(1+exp(mean_la*beta2(1:k)+U(RAND(j),1))+exp(mean_lb*beta2(k+1:2*k)+U(RAND(j),2)));
    
    PART(k+2:k+4,j)=exp(mean_hb*beta2(k+1:2*k)+U(RAND(j),2))...
        ./(1+exp(mean_ha*beta2(1:k)+U(RAND(j),1))+exp(mean_hb*beta2(k+1:2*k)+U(RAND(j),2)))...
        -exp(mean_lb*beta2(k+1:2*k)+U(RAND(j),2))...
        ./(1+exp(mean_la*beta2(1:k)+U(RAND(j),1))+exp(mean_lb*beta2(k+1:2*k)+U(RAND(j),2)));
end

part=PART*ones(draws,1)/draws;
